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Normal human heart rate shows complex fluctuations in time, which is natural, since heart rate 
is controlled by a large number of different feedback control loops. These unpredictable fluctuations 
have been shown to display fractal dynamics, long-term correlations, and 1/f noise. These char- 
acterizations are statistical and they have been widely studied and used, but much less is known 
about the detailed time evolution (dynamics) of the heart rate control mechanism. Here we show 
that a simple one-dimensional Langevin-type stochastic difference equation can accurately model 
the heart rate fluctuations in a time scale from minutes to hours. The model consists of a determin- 
istic nonlinear part and a stochastic part typical to Gaussian noise, and both parts can be directly 
determined from the measured heart rate data. Studies of 27 healthy subjects reveal that in most 
cases the deterministic part has a form typically seen in bistable systems: there are two stable fixed 
points and one unstable one. 

PACS numbers: 87.19.Hh, 02.50.Ey 



I. INTRODUCTION 
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Various methods and models have been used in at- 
tempts to characterize the dynamics of the heart rate 
control mechanism. For short time periods and under 
stationary conditions there are successful models of heart 
rate and blood pressure regulation Q, , but the charac- 
terization of long-term behavior has been a very difficult 
problem. Some models have been introduced in order 
to explain long-term fluctuations, but usually they can 
only describe well-controlled in vitro experiments, or the 
models depend on large number of parameters, which 
cannot be easily determined from experimental data 0. 
Furthermore, these models can predict only global sta- 
tistical features like scaling properties of power spectrum 
and correlations |j| and tell us very little about the de- 
tails of the time evolution. 

Many features can be extracted from long time se- 
ries of heart rate measurements, quantities like entropy 
measures 
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exponents [2TI |24| ) and symbolic dynamics complexity 
pa l2rj E3 , but these are all purely statistical charac- 
terizations and as such cannot provide us a mathemat- 
ical model of heart rate dynamics, not even a simple 
one. However, some of these statistical methods do char- 
acterize the complexity of the dynamics underlying the 
time series 28], or are directly related to their fractal or 
chaotic features. Mathematical analysis of many phys- 
iological rhythms, including long-term heart rate fluc- 
tuations, has revealed that they are generated by pro- 
cesses which must be nonlinear, since linear systems can 
not produce such a complex behaviour |29j. Nonlinear 
purely deterministic models can display chaotic dynam- 
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FIG. 1: Typical R-R interval time series recorded at night 



ics and generate apparently unpredictable oscillations, 
but in practice it has not yet been possible to extract 
such models from real noisy experimental data. It is also 
possible that the underlying system is stochastic, i.e., the 
time evolution of the system is subject to a noise source. 
[This kind of dynamical noise is different from measure- 
ment noise, which is mostly generated in the experimen- 
tal apparatus.] In any case, there is increasing evidence 
that noise, originated cither from the system itself or as 
a reflection of external influences, is actually an i nteg ral 
part of the dynamics of biological systems [3(|, |3ll |32] ■ 

A typical R-R interval recording is shown in Fig. ^ 
The time series is generated by recording a 24-hour elec- 
trocardiogram and detecting the R-peak from each heart- 
beat, the R-R interval is the time difference between two 
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FIG. 2: Schematic presentation of the method for analyzing stochastic time series and calculating the deterministic and 
stochastic part of the dynamics (Fig. A). Whenever the trajectory of the system passes near certain point x in the state phase, 
i.e. X(ti) ~ x, the future value X(ti + r) of the trajectory is recorded. The distribution of these values is fitted by a Gaussian 
function with the mean x + g(x) and deviation h{x), cf. equation This is repeated for all x-values. On the right are two 
typical examples of the distribution of future values, the initial x values are marked with dashed vertical bars and the fitted 
Gaussian curves with a thick line (Fig. B). 



consecutive R-peaks. In the upper panel of Fig.QJwe have 
the R-R interval time series for 6 hours. We can see sec- 
tions where the oscillations are rather regular but there 
are also abrupt changes. In the lower panel of Fig. ^ we 
have zoomed into a part of the time series, about 50 min- 
utes, and also on this time range we can see apparently 
random oscillations with rapid changes. 

It is well known that most short-time fluctuations of 
heart-rate are generated by respiration (periods typically 
in the couple of seconds range) and blood pressure reg- 
ulation (so called Meyer waves with periods of about 10 
seconds HH). In the following we are not interested in 
these fast rhythms (which can be analyzed quite well us- 
ing linear or semi-linear models) but rather in time scales 
from minutes to hours. We will show that in this time 
range the dynamics of the heart-rate fluctuations can be 
well described by a one-dimensional Langevin type dif- 
ference equation. This equation contains a deterministic 
part and additive Gaussian noise, and we have found that 
it works well when the delay parameter in the equation 
is in the range of 2-20 minutes. 



II. THE MODEL 

An important and wide class of dynamic systems can 
be described by the Langevin differential equation |34l 

M 



dX(t) 
dt 



g(X(t),t) + h(X(t),t)T(t). 



(1) 



Here X(t) represent the state of the system at time t, the 
function g gives the nonlinear deterministic change, and 



the last term h is the amplitude of the stochastic con- 
tribution and r(t) stands for uncorrelated white noise 
with vanishing mean. These kinds of stochastic differen- 
tial equations always need an interpretation rule for the 
noise term, normally one uses the Ito interpretation |36j| . 
In general the functions g and h could depend explicitly 
on time t. The equation (|TJ can be easily generalized to 
higher dimensions. We will now show that long-term be- 
haviour of heart rate can be modeled using a difference 
version of the Langevin equation |.31| 

X(t + r) = X(t) + g(X(t);r) + h(X(t)-r)T(t). (2) 

Here X(t) again represents the state of the system, in this 
case the R-R interval, at time £, and r is the time delay. If 
arbitrary small delays r are possible then one can take the 
limit t — ► and get the differential equation (JIJ [if the r 
dependence is given by g(X(t);T) « rg(X(t))], but in the 
present case it will turn out that there is a minimum r for 
which the model J2J seems to be valid. We assume that 
g and h do not have explicit time dependence but they 
may depend on the delay r. It is convenient to extract 
the term X(t) in the deterministic part, as is done in J2J, 
then a nonzero g(X(t);r) stands for changes in the state 
of the system. An essential feature of models of the above 
type is that for time evolution we only need to know the 
state at one given moment and not its evolution in the 
past, i.e., they are Markovian [34ll3^|. 

The computational problem is now to determine the 
functions g and h from measured time series and to ver- 
ify that the description using (0) is accurate. The prin- 
ciple of the method is very simple |38, 39] : at every time 
ti when the trajectory of the system meets an arbitrary 
but fixed point x in state space, we look at the future 
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state of the system at time f, + r. The set of these fu- 
ture values (for a chosen x and r) has a distribution in 
the state space and from this distribution we can deter- 
mine the deterministic part g(x) and the stochastic part 
h(x), see Fig.[2\. In practice we first divide the range of 
the dynamical variable X into equal boxes. By scanning 
the whole measured time series we check when X is in- 
side a given box x, i.e., \X(ti) — x\ < Ax , where x is 
the middle value of the box and Ax is the half width of 
the box. When X is found on the box, we look at the 
future value of the variable, X(i, + t), where r is the 
fixed delay parameter. Since the trajectory of the sys- 
tem passes each box several times, we can calculate the 
distribution of the future values X(U + t) for each box 
x. If we assume that the noise is Gaussian, we can fit a 
Gaussian function on each distribution, and as a result 
we get the mean and the deviation parameters for each 
x; the mean of this distribution is equal to x + g{x) and 
the deviation is equal to h(x) |4(J, HJ. A typical case 
is given in Fig. [23, and it shows that the distribution is 
actually very well described by Gaussian noise (the cor- 
relat ion is better t han 0.95; the correlation is calculated 
as yj\ — S res /Stot, where S res is the sum of the squared 
residuals and Stot is the variance). From the given data 
we can in this way determine the functions g(X) and 
h(X) needed in the stochastic model (J5J. It should be 
noted that we can calculate only the absolute value of 
h(X) since the deviation parameter found from the fit- 
ted Gaussian function is in squared form. 

In our analysis we have used R-R interval time series 
of 22-24 hours, corresponding to 80.000-100.000 data 
points. Our data is actually interval data, i.e., it consists 
of a sequence of R-R interval values. It is then convenient 
to count the delay in our analysis in terms of heart-beats 
rather than seconds, i.e., we have not used cumulative 
time as time variable but the beat index. However, since 
the R-R interval values vary a lot within the used delay 
range, the beat index actually gives a delay as if com- 
puted with the average beat-rate. We have tested both 
methods and found only minor differences between them 
(in the details of the functions g and h) . We will show 
later that the functional forms of g and h are quite in- 
sensitive on the time delay, and since this holds for both 
methods we will use the more convenient beat index. 



III. RESULTS 

A. A typical case 

In Fig. |2| we have presented results obtained for a par- 
ticular case using the method described earlier. The 
value of the delay parameter r was 500 beats, and the 
number of boxes used to construct local distributions was 
150. Distributions were fitted using Gaussian function. 
The g{X) function, the deterministic part of the system, 
is displayed on the left panel in Fig 3. It has a very 
clear and simple functional form (between the vertical 



lines) which is typical for systems exhibiting bistable be- 
haviour 0, ^3 • The function crosses the zero line three 
times, these crossings are the fixed points of the system. 
The fixed points marked with arrows are stable: with- 
out any noise term these points attract all nearby states, 
because the control function g(X) is locally decreasing. 
The middle fixed point is repulsive. Due to the stochastic 
part the system has a tendency to jump between the sta- 
ble points if the amplitude of the noise is high enough. 
Far away from the stable points g{X) increases or de- 
creases strongly forcing the system rapidly back to os- 
cillate around the stable points. The amplitude of the 
stochastic part of the system, function h(X), is almost 
constant except between the stable points where it has 
a clear maximum (the middle panel in Fig. EJ) . One in- 
terpretation is that the system has a larger inherent free- 
dom to oscillate randomly when the trajectory is between 
the stable points but outside this range the character of 
the system is more deterministic. From the physiological 
point of view this kind of dynamics can be useful since 
it lets the R-R interval to wander most of the time but 
prevents it from escaping too far away from the normal 
range. On the right panel in Fig. |3| we have shown the 
correlation coefficient of each local distribution. Most of 
the time the correlation is remarkably high, about 0.85- 
0.95, but near the largest and smallest X values there 
are only rather few data points and therefore the corre- 
sponding distributions do not have clear Gaussian shape 
resulting with lower correlation. The high average cor- 
relation value is a clear indication that the noise in this 
system is really Gaussian type. We have used the value 
of 0.8 as a threshold level, and the corresponding range 
is marked with the vertical lines in Fig. |3J 

What is remarkable in this description is that the func- 
tional forms of g(X) and h(X) are fairly independent of 
the delay parameter r in a rather extensive delay range, 
typically 100-1000 beats (corresponding to 2-20 min- 
utes). In Fig. 01 we have plotted the functions g(X) and 
h(X) for a range of r values. The ^-function is practi- 
cally r independent, except for the shortest R-R intervals, 
where some cumulative effects show up. The ft-function 
seems to grow very slowly as r increases. For still smaller 
delay values g(X) is more flat and h(X) is more scat- 
tered, and for longer delays g(X) is typically a straight 
line and h(X) is constant. Behavior at these extremes 
can be easily understood by recalling that when the time 
scale is small, the heart rate system is clearly multidi- 
mensional depending directly on blood pressure, respira- 
tion and other rapidly changing physiological variables 
and our 1-dimensional description is no longer valid. On 
the other hand, if the delay parameter is very large, we 
cannot reconstruct the local dynamics in terms of local 
distributions, we just get the global distribution that is 
independent of dynamics and no longer Gaussian. In the 
right panel of Fig.0]we have given the values of the g(X) 
function at X = 780 ms (marked with a vertical dashed 
line in the left panel) computed with delays of 5-10240 
beats. We can see a plateau in the delay range of 100- 
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FIG. 3: Typical results derived from R-R interval time series using time delay r = 500. We have shown the deterministic part 
g{X) (the left panel), the stochastic part h(X) (the middle panel) and the correlation coefficient of the distribution (the right 
panel) as a function of the dynamical variable X. The range corresponding to the correlation threshold level of 0.8 is marked 
with the vertical lines. 
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FIG. 4: Examples of the deterministic part g(X) (the left panel) and the stochastic part h(X) (middle panel) calculated with 
various values of the delay parameter r : 40 (thick line), 80 (thin line), 160 (thick line), 320 (thin line) and 640 (thick line). 
The values of the g(X) function at X = 780 ms (marked with vertical dashed line in the left panel) are plotted as a function 
of the delay in the right panel, there is a plateau around a delay of 100-1000 beats. 



1000 beats which means that the g(X) curves for these 
delays are bundled. In principle the curves for a delay 
of 2t should be obtainable by iterating J2J with delay r. 
Direct numerical calculations of joint probabilities using 
experimentally determined g(X) (within 100-1000 beats 
delay range) indicate that g(X) and h(X) do not change 
significantly in one iteration, mostly because in our case 
the Gaussian distribution is not so narrow. In general it- 
erations tend to sharpen the bends in g(X) and this fea- 
ture is indeed visible in Fig. The small r-dependency 
of g and h in the range of short R-R intervals can then 
be interpreted either as the expected result from repeated 
iterations, or as a sign of higher order dynamics: possibly 
the heart rate regulation system is more complex when 



the system must readjust at a fast heart rate. 



B. Variation between subjects 

In order to find whether different subjects have are 
any common features in the deterministic and stochas- 
tic parts g(X) and h(X) we analyzed the data from 27 
healthy subjects of various age and gender [18 cases from 
PhysioBank 0] and 9 cases from Kuopio University Hos- 
pital]. Analyses were done using the same parameter 
values as in Fig. [3J The deterministic part, the g(X) 
function, is displayed in Fig. [S] for a set of 9 typical cases. 
The most common form for this function is the bistable 
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FIG. 5: Typical deterministic functions g(X) derived from different subjects. The cases 1-5 represent simple bistable situation, 
the cases 6 and 7 have three stable points, the case 8 is multistable, and the case 9 has only a single stable fixed point. 



type, already shown in Fig. EJ where the g{X) function 
has three zeroes, and 60% of all cases can be classified to 
this group (cases 1-5 in Fig. EJ ■ The next most common 
group, 25% of all cases, has a g(X) function with 5 ze- 
roes, a kind of double pitchfork system (cases 6 and 7 in 
Fig. EJ). We also found 3 cases where the g{X) function 
seems to have even more zeroes (case 8 in Fig. EJ . Only 
very few cases could not be clearly classified as bi- or 
multistable. In these cases it can be difficult to interpret 
the results. It is possible that the dynamical variable did 
not explore the whole state phase, and therefore we can 
see only part of the g(X) function; for example case 9 in 
Fig. El where the system has only one stable fixed point 
and no unstable points at all, can be an example of this. 



The stochastic parts (function h(X)) are fairly similar: 
they are almost constant except that in all cases there are 
maxima on the R-R interval ranges between the stable 
fixed points of the deterministic part, as in the example 
in Fig. E| 

The description given by equation J2J contains both a 
deterministic and a stochastic component. It is an im- 
portant to realize that the stochastic part is not a small 
perturbation but in fact forms an essential part of the de- 
scription, furthermore it is 10-20 times higher than the 
measurement noise (uncertainty in detecting the position 
of the R-peak), which is typically only 2-5 ms. One way 
to compare the deterministic and stochastic components 
is to note that the size of the bend in the g(X) function is 
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FIG. 6: The deterministic parts g(X) (left panel) and stochastic parts h(X ) (right panel) computed from the R-R interval time 
series recorded from the same subject on different days. The data from the first recording is marked with solid dots and data 
from the second recording, 4 days later, with open dots. 



of the order of 30 to 50 ms, while the average size of the 
h(X) function is about 70 to 110 ms, as can be seen in 
Fig- El [The extraction of small details in the g(X) func- 
tion under such noise is of course possible only because 
the noise is so cleanly Gaussian.] On the other hand, 
the distance between the stable fixed points in the g(X) 
function is of the order of 50 to 250 ms, and therefore the 
probability that the systems jumps between stable points 
is not extremely high, but nevertheless possible. It is also 
possible that external factors drive the system from one 
stable point to another, since during night-time the mean 
R-R interval is typically longer than during day-time [al- 
though the R-R interval can abruptly jump to the faster 
rate also during the night, as can be seen on the lower 
panel in Fig. P . 



C. Same subject at different times 

If the model Q were to describe true heart rate dy- 
namics, the functions g(X) and h(X) should have some 
constant features specific for each subject. In order to 
look at this aspect we made two recordings from the same 
subject within 4 days, the results are shown in Fig.El In 
general the deterministic and stochastic parts from differ- 
ent recordings are remarkably similar both having clear 
bistable character. In the R-R interval range of 500-800 
ms the results are almost identical and the only difference 
seems to be a scaling towards the shorter R-R intervals 
in the 800-1100 ms range of the second recording. In the 
first recording the mean value of R-R interval calculated 
over the 24 hour period was 781 ms and in the second one 
726 ms. Therefore in the second recording the shortest 
R-R intervals are significantly more frequent and this can 
affect on the analysis results. These deviations could also 
reflect true changes on the underlying control system: it 
is well know that there are daily variations on functions 
of the autonomic nervous system. 



D. Surrogate analysis 

As a further validity check we also performed surrogate 
analysis |4J, |45j in order to eliminate the possibility that 
the results are generated just from peculiar distribution 
of the R-R intervals imitating real dynamics. For this 
purpose the data was shuffled by dividing it into sections 
of equal size which were then repositioned randomly. As 
a result we get a new time series where the dynamical 
structure has been partially destroyed depending on the 
section size. Results of this surrogate analysis are shown 
in Fig. The top panels display the deterministic g(X) 
and stochastic h(X) parts of the system and the corre- 
lation coefficient without any data shuffling (row A in 
Fig. UJ. On the next row (row B in Fig.0) we have used 
sections of 800 data points for shuffling. There are only 
small changes in the deterministic part, but the correla- 
tion has decreased noticeably. When the section size is 
400 (row C in Fig. [7J we can no longer see the bistable 
character in the deterministic part, the stochastic part is 
flat with higher mean level, and the average level of the 
correlation coefficient has dropped well below our thresh- 
old value 0.8. With still smaller section sizes the results 
do not change any further. In this analysis we have used 
the same delay of 500 data points as previously and when 
the section size used in the shuffling process is less than 
this delay all dynamical properties disappear, as expected 
in the case of true time evolution. Therefore we conclude 
that our results are derived from the dynamical proper- 
ties of the heart beat data, and not from their overall 
statistical characteristics. 



IV. CONCLUSION 

Our results indicate that the human heart-rate con- 
trol dynamics can be accurately modeled with the 1- 
dimensional stochastic difference equation J2J , where the 



7 



h(X) 



Correlation 



B 




1050 1400 700 

X[ms] 



1050 

X[ms] 



1050 

X[ms] 



FIG. 7: The deterministic part g(X) (left column), stochastic part h(X) (middle column) and correlation coefficient (right 
column) for the original data (row A) and for two surrogate versions (B and C). For surrogate data the original data has been 
shuffled using section sizes of 800 (B) and 400 (C) data points. 



time delay parameter is within 2-20 minutes. Stochas- 
ticity is an integral part of the dynamics, and in this 
delay range the effects of other variables are either em- 
bedded into the stochastic part of the system or averaged 
over time with no net effect. It is remarkable that the 
form of the control function g{X) is similar from case 
to case. Their typically bistable character is also well 
justified on common physiological grounds. From this 
initial study we cannot yet identify what kind of dynam- 
ical structure is typical for healthy subjects (although our 
results already indicate that a simple bistable system is 
most common feature) and therefore the model cannot 
yet be used directly for clinical work, for that purpose 
one needs extensive demographic studies. We can nev- 
ertheless speculate that the form of the control function 
g(X) should tell us something about the health of the 



subject. Also, some of the current knowledge based on 
statistical measures of heart rate time series can prob- 
ably be explained within the framework of our model. 
Another interesting observation is the importance of the 
stochastic part, it could be the result of integrating the 
effects of a more detailed control mechanism over time, 
but it could also reflect some truly stochastic internal and 
external influences. 
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